Exact evolution of time-reversible symplectic integrators 
and their phase error for the harmonic oscillator 
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The evolution of any factorized time-reversible symplectic integrators, when applied to the har- 
monic oscillator, can be exactly solved in a closed form. The resulting modified Hamiltonians 
demonstrate the convergence of the Lie series expansions. They are also less distorted than modified 
Hamiltonian of non-reversible algorithms. The analytical form for the modified angular frequency 
can be used to assess the phase error of any time-reversible algorithm. 
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I. INTRODUCTION 

Symplectic integratorsii are methods of choice for 
solving diverse physical problems ranging from celestial 
mechanics^, molecular dynamics'^ , to accelerator physics^ 
and lattice gauge theory?:. In contrast to other numerical 
methods for solving Hamiltonian dynamics, symplectic 
integrators evolves the system according to a modified 
Hamiltonian^-^ that can be made arbitrary close to the 
original Hamiltonian. For this reason, despite the fact 
that symplectic integrators can preserve all Poincare in- 
variants, they can never exactly conserve energjii. Oth- 
erwise, an integrator's modified Hamiltonian would have 
coincided with the original Hamiltonian, and the integra- 
tor's evolution would have been exact. 

Symplectic integrators can be derived most easily by 
the method of factorization/composition^-S, where the al- 
gorithm is composed of many elementary ones with vary- 
ing coefficients. In this case, one can systematically com- 
pute the modified Hamiltonian order-by-order via the Lie 
series^ approach by using the Baker-Campbell-Hausdorff 
(BCH) expansionSi^. While this approach can demon- 
strate the existence of the modified Hamiltonian pertur- 
batively, and show that it can be made arbitrary close to 
the original Hamiltonian, one has no sense of the modified 
Hamiltonian's global structure nor of the series' conver- 
gence. 

In this work, we show that the evolution of any factor- 
ized, time-reversible integrators, can be exactly solved 
for the harmonic oscillator. The algorithm's evolution 
remained that of a harmonic oscillator but with modified 
mass and spring constant. Thus the model illustrates 
very simply, how reversible symplectic algorithms pre- 
serve invariant tori of a dynamical system. Moreover, 
the closed form expression for the modified Hamiltonian, 
when expanded, matches the Lie series order by order, 
thus demonstrating the latter's convergence. Finally, we 
show that the analytical modified oscillator frequency can 
be used to benchmark any time-reversible algorithm. 

In the following, we will briefly summarize key features 
of the Lie-Poisson construction of symplectic integrators, 
the Lie series expansion for the modified Hamiltonian, 
the phase-space matrix approach, and comparative re- 
sults for a number of fourth order reversible algorithms. 



II. SYMPLECTIC INTEGRATORS 

The evolution of any dynamical variable W{qi,pi), is 
given by the Poisson bracket, and therefore by the cor- 
responding Lie operator H associated with the Hamilto- 
nian fimction H{qi,pi), i.e. 



dW 
dt 



{W,H} 
/dH d 



dH d 



V dpi dqi dqi dpi 



W = HW. 



(1) 



For any dynamical function Q, we can similarly define its 
associated Lie operator Q via the Poisson bracket 

QW^{W,Q}. (2) 

The operator equation can be formally solved via 



Wit) 



W{0). 



(3) 



Symplectic algorithms are derived by approximating the 
evolution operator e*^ for a short time in a product form. 
For Hamiltonian function of the usual separable form, 

i7(q,p) =T(p) + y(q), with T{p) = ^p,p,, (4) 

the Hamiltonian operator is also separable, H = T + V, 
with differential operators T and V given by 



f = {.,r} 



dT d 
dpi dq.i 



Pi 



d_ 

dqi 



dqi dpt 



d 

OPi 



(5) 



(6) 



The Lie transforms^ e^"^ and e^^ , are then displacement 
operators which shift qi and pi forward in time via 



and 



(7) 



Thus, if e^^ can be factorized into products of Lie trans- 
forms e'^-^ and e^^ , then each factorization gives rise to an 
integrator for evolving the system forward in time. Ex- 
isting literature on symplectic algorithms are concerned 
with decomposing e'^^ to high orders in the product form 



N 

n 



(8) 
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with a set of well chosen factorization coefficients {U, Vi}. 
In most cases, we will consider only (left-right) symmetric 
factorization schemes such that either ti = and Vi — 
VN^i+i, U+1 = tN-i+1, or vn = and = VN-i, U = 
tAT-i+i. In either cases, the algorithm is exactly time- 
reversible, and the energy error terms can only be an 
even function of e. Such a symmetric factorizations is 
then at least second order. 



III. LIE SERIES EXPANSION 

A well known symmetric algorithm is the second order 
Stormer/Verlet (SV) algorithm, 

Tsv{e) = e^^^e'^e^'^ = e'"^ , (9) 

where Ha is the approximate Hamiltonian operator of 
the form 



HA=f + V + e^ (e^^^ [f^V] + e 



VTV 



-\-e^ 



[T{fvf 



.[V{TVf 



., (10) 



where e^^j,, e^r^,^^^ etc., are coefficients specific to a 
particular algorithm. We have used the condensed com- 
mutator notation [f'^V] = [f,[f,V]]. For the case of 
algorithm SV, the error coefficients are 



_ 1 
_ 1 

^TTTTV 720' 

_ 1 

"ttvtv "120' 



1 

24' 

_ 1 

~ 360 
1 

480 



-7^- (11) 



These are computed by repeated applications of the 
Baker-Campbell-Hausdorff (BCH) formula for combining 
exponentials of operators. Thus the algorithm evolves 
the system according to Ha rather than the true Hamil- 
tonian H. Knowing Ha then allows us to determine the 
actual Hamiltonian function Ha which governs the algo- 
rithm's evolution. From the fundamental definition of Lie 
operator (PI, one can convert operators back to functions 
via [T,V] ^ {V,T} = -{T,V}. For symmetric factor- 
ization with only even order commutators, the Hamilto- 
nian function can be obtained simply by replacing con- 
densed commutators with condensed Poisson brackets: 
[T^y] ^ {T'^V} = {T,{r,y}}, etc.. For the harmonic 
oscillator, 



(12) 



the non- vanishing Poisson brackets are: 

{VTV} = 



2 2 
u) p , 



{TTV} 
{T{TVf} ^ -2cjV, 



4 2 
-LU q , 



There is thus a clear separation between the error coeffi- 
cients, which are characteristics of the algorithm and the 
Poisson brackets, which are properties of the underlying 
Hamiltonian. If the SV algorithm is applied to the har- 
monic oscillator, then the resulting modified Hamiltonian 
function is 



HA{q,p) 



1 



1 



■P 



1 * 2 
k q 



2 m* 2 
with effective mass and spring constant 



1 

TO* 


= 1 + 


6 


1 4 4 , 

'30^" 


k* 


= 1 - 


loo 

12^" 


1 44, 
-120^" + 



and approximate angular frequency loa = 
algorithm therefore evolves the system according to 

sin(uj At) 



(14) 

(15) 
(16) 

The 



f q{t) \ _ cos{uJAt) 
[pit) 




(17) 



— \/m*k* sm(LLjAt) cos{u!At) 
By contrast, a first order, non-reversible factorization 

e-V^ = e^^-, (18) 
with Ha directly given by the BCH formula. 



Ha - H+^e[fv] + ^e[ffv] - ^e[Vfv] - 
produces the following expansion. 



(19) 



Ha = ^[p^ +uj^q^ +eLu^pq 



30 



4 4 
e uj 



(20) 



{ViTVy} = 2Lu'>q\ (13) 



Note that the overall factor is the same as (|15ll . The 
odd order terms distort the harmonic oscillator funda- 
mentally, producing a distinctive, 45°-tilted phase-space 
ellipse, commonly seen in the standard map^, or when 
solving the pendulum^. This is an artifact of not being 
time-reversible, a property foreign to the original dynam- 
ics. Time-reversible algorithms with modified Hamilto- 
nian such as H14|l showed no such distortioniS. 



IV. MATRIX METHOD 

The above expansion method for computing the modi- 
fied Hamiltonian, while general, is severely limited to low 
orders. At higher order, the number of commutator and 
Poisson brackets proliferates greatly due to the complex- 
icity of the BCH expansion, it is then difficult to assess 
the convergence of series (fT5|) . (fTB|l and ((201 ■ 

To compute Ha exactly, we completely abandon the 
BCH-based approach. For the 1-D harmonic oscillator. 
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let's denote the phase-space vector at the ith iteration of 
the algorithm as 




For the SV algorithm, let T and V denote the effect 
of Lie operators e*^"^ and e^^^ acting on these vectors. 
From Q, it is easy to see that they are upper and lower 
triangular matrices given by 

Hence, the algorithm corresponds to the product matrix, 
M = VTV = ( ^'Y2 2^1 f 2 2 V (23) 

Matrix representations have been widely used to study 
the stabilitjiiiii^ and convergenc o^i^'^ of low order sym- 
metric integrators. Our interest here is not to solve the 
harmonic oscillator per se, but to use it as a mean of ob- 
taining the converged value of the Lie series. In contrast 
to previous use of the matrix method to study the har- 
monic oscillatori2iiiii24i^, we are not interested how each 
algorithm's approximate solution approaches the exact 
solution; we are interested only in the exact form of the 
approximate solution itself. Below, we also emphasize the 
distinction between time-reversible and time-irreversible 
algorithms. As illustrated in H22|l . for quadratic Hamil- 

tonians, each elemental operator e''^"^ and e"'^^' can in 
general be represented as upper or lower triangular ma- 
trices, 

where the constant matrix elements (7i{s) and /ii(e) are 
both positive as e — > but are odd functions of e: 

(Ti(-e) = -cTiie), f^i{-e) = -fj.i{e). (25) 

Both Ti and can be further decomposed as 

T.(£)-(^J ?)+^'(^)(!! l)=i + Me)t, (26) 

V.(£)=(^J ?)-A^.(e)(? q) =1-A^.(£)v, (27) 

where t and v are nilpotcnt matrices, 

= 0, = 0, (28) 

with property tv + vt = 1. The nilpotence of t and v 
makes it obvious that both and are time-reversible, 

T,(-£)T,;(£) - 1, V,(-e)V,(e) = 1. (29) 



even in s and off-diagonal elements odd in e. In fact such 
a general matrix 

Gie)^gie')l + r{e)t-^ie)v=(^^^ ^) (30) 

with t(— e) = —t{s) and e) = —v{e), is indeed time- 
reversible, 

G(-£)G(£) = g^l-{Tt-iyvf 

= (5^ + iyr)l= (dctG)l, (31) 

provided that (det G) = 1. Since this is true of any 
symplectic integrator, we conclude that any symplectic 
integrator of the form (|30|l . is time-reversible. The gen- 
eralization of (|23|) , corresponding to the symmetric form 
of can be now written as 

M(£)=nT.(£)V,(e). (32) 

i 

The advantage of the matrix representation is clear. Each 
matrix product can be computed exactly, whereas the 
corresponding product of exponential operators would 
have required an infinite BCH expansion. For symmet- 
ric factorization, M is time-reversible. We now prove by 
induction that M is indeed of the form (|30|) . Let G be 
of the form (|30|l . then 

G' = T,GT„ 

= gT^ -I- rT.tT, - z^T.vT,, 

= {g - va,)t+ {t + 2ga., - iy(7f)t - vv, (33) 

which is again of the form (|30|l . Similarly for V^GVi, 

G' = ,gV2 + rV.tV, - i.V,vV, 

= {g-TpL,)l + Tt-{,^ + 2gfi,^Tfi'^)v. (34) 

Since the initial G is either or V^, which is of the 
form H3U|I . the proof is complete. If the algorithm is not 
time-reversible, then the diagonal elements will not be 
equal and M will have the more general form 

M.(4, ;;). (35) 

We will consider this more general matrix and recover 
the time-reversible case at the end by setting h = g. 

For the harmonic oscillator, because the matrix M is 
constant for all time steps, we have 

r„ - M^ro. (36) 

To evaluate analytically, we diagonalize M via 

7W = S-iMS. (37) 

For VT > (g — /i)^/4, the eigenvalues are unitary. 



The time- reversibility of Ti(e) and Vi(e) can also be 
traced to the fact that they have equal diagonal elements 



(38) 
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where 



^ = ^/vT - [g - hf/A, 
'1 



cos 



■TrM 



9- 



2 J V 2 

The diagonal matrix and the transform are therefore 

-ie Q \ / (h-g)/2+ii {h-g)/2-i, 



(39) 
(40) 



M 







1 



1 



The A^th power of M can now be computed via 
where we have substituted N — t/e, 



(41) 



A = -{ie)s 

£ 



-1 
1 



2 



70* lit 

-72 t -7o t 



yielding finally 

9 [g-h) 

70 - 



71 



y r 



72 



(42) 
(43) 

(44) 



£ 2^ 

The exponentiation of A then gives, 

M^ = R+E, (45) 
where R is a pure rotation and S, a translation: 
cos(f) fsin(f)- 
_^sin(^) cos(^^) 



R 



{g-h) . (9t 
S = ^^sm(- 



1 

-1 



(46) 
(47) 



For time-reversible algorithms, /i = g, S = and ^ = 
y/vT. By comparing H46(l to (|17|l one can quickly identify 



UJA 



UJA 



e_ i_ 

e ' V V T 



k* — LOA 



(48) 



and the modified Hamiltonian conserved by the algo- 
rithm. 



(49) 



The analytical form of the algorithm's solution (|46|l is 
interesting in that the discrete character of the algorithm 
has disappeared! The solution is a continuous function 
of i, not necessarily a multiple of e. 

For non- reversible algorithm such as H18|) . S 7^ 0. This 
operator translates q ^ q + s and p ^ p — s peri- 
odically, thus causing an elongation of the phase-space 
ellipse along the minor diagonal at -45°. This is the 
phase-space distortion alluded to earlier and observed 
in the standard map^. While non-reversible algorithms 
may be of interest mathematically, from the perspective 
of physics, there is nothing to be gain in solving time- 
reversible dynamics with non-reversible algorithms. 



V. MODIFIED HAMILTONIANS AND THE 
PHASE ERROR 

To illustrate how this formalism works, for the SV al- 
goritm, we have from (|23|) . 



1 1 2 2 



£, = £w2(l - -£^^2). (50) 



From these we can determine in closed forms 
UJA = icos"^(l - ^£^uj^), 

m* UJ 4 
k* = UJAUJ(1 - ^£^^^2)^/2. 



(51) 
(52) 
(53) 



We can now expand these physical quantities to any or- 
der, but for the lack of space, we will stop at the 6th: 



+ •••, (54) 
(55) 
(56) 



UJA 


- 1 + 




3£4w4 


5£^w6 


UJ 




24 


^ 640 


7168 


1 


= 1 + 


e^uj^ 


4 4 
e uj^ 


e^uj^ 


m* 


6 


+ 30 ^ 


140 ^ 


k* 


= 1- 


£2.^2 


4 4 
e UJ* 


e^uj^ 


Z7 


12 


120 


840 ^ 



The first three terms in each case are in agreement with 
(I15|l and H16|l . By comparing (|54(l to 1)51(1 one can im- 
mediately conclude that for < £0; < 2, the Lie series 
converges. Thus the Lie series converges up to the in- 
stability point. For another example, in the case of the 
fourth order Forest-Ruth integrator^'*, (We use factoriza- 
tion coefficients as derived by Creutz and Gocksch^ and 
Yoshida^^) 



9fr 



1 - 



288 



(144 - Ue'^uj^ - (6 + 5 v/2 + 4^)£ 



^e^uj^ 



6 72V22 
^(25 + 20^4- 16 V^)£<'w6 



1728 



= euj 



6 144 



From these we can determine uja, l/m-¥ and k* via H48|l . 
They can then be expanded to any order, and for lack of 
space, we will stop at the 6th: 



UJA 

UJ 



(32 + 25v^+20v^)£4t. 
1440 

{S% + ^Q^/2 + bQ^/2?)£^UJ^ 
24192 

J_ _ ^ (6 + 5^+5^)£^cj^ 

m* ^ " 720 

(71-h56^ + 42v^)£6w6 

12096 



5 



k* _ ^ (26 + 20^+15v^)e'^w4 



720 



(80 + 63^ + 49sy22)e6cj6 
6048 

The modified Hamiltonian now deviates from the original 
Hamihonian beginning at the fourth order. 

In solving for periodic motion, the most serious error is 
not the energy, which is well conserved by symplectic in- 
tegrators after each period. The most serious error is the 
phase (or angle) error, which accumulates after every pe- 
riod. The phase error thus grow linearly with integration 
time, even for symplectic integrators. The phase error 
is simply related to the fractional deviation of the the 
modified angular frequency from the original frequency: 



A(p = (lua — ti')T : 



UJ 



I) 



(57) 



For an nth order algorithm, its phase error is essentially 
given by the coefficient c„, 



IMA 
UJ 



1 



1 



euj 



+ 0(e"+^) (58) 



For the SV algorithm, we have from (|54|l . C2 — 1/24. 
Below we list the C4 coefficient for some recent fourth 
order algorithms. 



Ci{M) 

C4(C) 



(32 + 25v^+20v^) 

1440 
-0.0661431 

-2956612 -M24595V47T 

2797262640 
-0.0000902971 

-0.0000133432 
^ = 0.000130208 



7680 



(59) 



FR is the Forest-Ruth algorithm which requires only 3 
force evaluations. M is McLachan's algorithm^^. with 4 
force evaluations. BM is Blanes and Moan's algorithmii 
with 6 force evaluations. C is Chin's all positive time 
steps, forward algorithm.^^ C which uses 3 force and 1 
force gradient evaluation. Note that coefficient is posi- 
tive only for the forward algorithm. Since for the same 
number of force evaluations, one can update the FR al- 
gorithm twice at half the step size as compared to the 
BM algorithm, one should reduce FR's error coefficient 
by a factor 2* = 16 before comparing it to BM's error 
coefficient. Or, one should multiplying BM's error by 
(6/3)'' before comparing to that of FR's. Thus for an 
equal computational effort comparison, we should multi- 
ply each algorithm's coefficient by (n/3)'' where n is the 
number of force evaluation and normalize it by dividing 
by FR's value. This normalized coefficient c* for each al- 
gorithm is then as follow, (For this comparison, wc count 



each force gradient as roughly equivalent to one force 
evaluation. For most forces in celestial mechanics and 
molecular dynamics, the force gradient can be evaluate 
without too much effort.) 



cliFR) 
cl{M) 
cl{BM) 



-1.0000 
-0.0043 
-0.0032 
0.0062 



(60) 



Since the dynamics of the harmonic oscillator is very 
simple, its phase error is only a necessary but not a 
sufficient criterion for gauging an algorithm efficiency. 
While it is well known that FR has relatively large er- 
ror among fourth order algorithms, it is interesting that 
this simple analytical calculation can indeed distinguish 
the latter three algorithms as significantly better than 
FR. (In the case of the forward algorithm, we have 
deliberately chosen algorithm C to do the comparison. 
Within the parametrized family of fourth order forward 
algorithmsiSiSS, one can actually make C4 vanish. Thus 
a general fourth order forward algorithm can solve the 
harmonic oscillator to six order, but such a comparison 
would not have been fair without allowing other algo- 
rithms to be fine-tuned.) 



VI. CONCLUDING SUMMARY 

In this work, we have shown that the evolution of 
any time-reversible symplectic algorithm can be solved 
in closed form. We proved that any factorized, time re- 
versible symplectic algorithm must have a phase-space 
matrix of the form H30|) . with equal diagonal elements. 
Based on the exact solution, we have further shown that 

1) the Lie series expansion for the modified Hamiltonian 
converges in the case of the harmonic oscillator. The 
global structure of the modified Hamiltonian is known. 

2) The phase-space structure of the harmonic oscillator 
is closely preserved by time-reversible algorithms but is 
distorted by non-reversible integrators. 3) The analyti- 
cal form of the phase error can be used to broadly assess 
the effectiveness of any reversible algorithm without do- 
ing any numerical calculation. We verified this claim by 
computing the phase error coefficients for four fourth- 
order algorithms. 
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